Just to say that I'm not sure if this falls under fix or feature. Clearly the particular use case of this update is fixing unintended behaviour when solvating, but I think it will be included with a more general update to the BioSimSpace.Solvent code that builds of OpenBioSim/biosimspace#149 Crystal waters seem to be a recurring issue at the moment, so I want to add an example page showing how to deal with them. I'm happy with whatever you decide for this particular PR.

Everything works as expected in BioSimSpace. Phew!

Actually, it's not working properly because the Sire GroTop parser doesn't seem to be able to handle multiple types of water molecule in the same file, so they either all get converted to the default template SOL, or that of the crystal water XTL, depending on the order of the molecules. Presumably this is because it does some matching for the parameters to avoid duplication.

This will need a rethink.

I guess this can be solved by changing one of the parameters for the molecule that won't affect the calculation by genion, e.g. something like the mass. Will try that to see if it works.

Nope, that doesn't work either, since it seems to get the mass from the element.

Hmmm, nothing seems to work, and I get a different topology file depending how the different water molecules (crystal and normal) are arranged in my system. Perhaps we'll need to use entirely different water templates for the crystal waters. If so, it's probably easier for me to just use a different approach for handling the genion issue, e.g. re-running until no crystal waters were removed.

This has also made me realise that the request in OpenBioSim/biosimspace#148, i.e. not changing the existing water topology on write on solvation, is not possible with the current GroTop parser, i.e. if there are other water molecules in the system then the molecule types will always be de-duplicated. (This is what you'd want for the purposes of running anything with GROMACS.) You should always be able to reference your original waters by index anyway.

The key line that detects water molecules in a GroTop file is here:

In theory, we could add an extra property, e.g. is_crystallographic and then change this line to read

auto waters = system.search("water and not molecule property is_crystallographic")

Then the crystallographic waters will be treated as separate objects, and won't be combined with all of the other water molecules.

The mols["water"] search string is very powerful, and will 100% find anything that looks like a water molecule. It does this by looking for molecules that contain one oxygen, two hydrogens and any number of ghost atoms. It is not based on mass or anything like that.

Another approach could be to change the elements or do something that means that mols['waters"] cannot find the water.

Or, to go deeper, we could add a crystallographic_water search term, e.g.

xtal_waters = mols["crystallographic_waters"]

and then we could use

auto waters = system.search("water and not crystallographic_water")

What approach do you think would be best?

(or, we could add a way of marking a water molecule as something that should not match the mols["water"] search term. This could be a is_not_searchable_water property?)

Aha, I didn't realise that you were now using the water search functionality. Yes, it should be pretty robust. I wrote it because I needed to be able to detect and swap water topologies in BioSimSpace, and this is obviously painful in Python for any reasonable number of water molecules.

In this case the crystal waters will have a very specific formatting since this is something that is only implemented as a fix behind the scenes. It should be simple enough to combine the water search with a residue name match to extract them, e.g.

xtal_waters = mols["water and resname XTL"]

I'll have a play with this tomorrow.

Thanks for the pointers!

That said, having the ability to flag something as not searchable might be quite useful, e.g. by adding a specific property to a molecule. This would help in the situation where we do want to preserve some user naming for a topology, but don't know what that naming would be, i.e. don't have an easy way to search for it.

Great, that works perfectly, i.e. adding excluding crystal waters from the search, then adding a separate search to define a new molecule type for them. I'll see how easy it is to add in a check for an is_non_searchable property, so that the user can keep custom water naming if desired.

I somehow hadn't noticed that the water searching that I wrote with IDEngineWater has now moved into functionality within SireMol itself, e.g. an is_water function. That's much more useful.

Okay, this all appears to be working. The BioSimSpace feature branch that uses these changes is here, along with an associated unit test.

Thanks, I agree. While this seems to work and is entirely invisible to the user, it would be good to have some time to test it in the wild and refine as necessary.